#include <stdlib.h>
#include <limits.h>
#include <fstream>
#include <vector>
#include <iostream>
using namespace std;

inline double rnd();
long double   binomialDistribution(const int k, const int n, const double p);

int main(int argc, char *argv[])
{
  const int      len = 40; // число "лунок" - 1
  const double   p   = 0.7;
  vector<double> probability(len);
  vector<double> cumulativeProbability(len - 1, 0.0);
  vector<long>   numberOfHits(len, 0);

  // инициализация
  for (unsigned i = 0; i < probability.size(); ++i) {
    probability[i] = binomialDistribution( i, probability.size(), p );
  }

  cumulativeProbability[0] = 0.0;
  for (unsigned i = 1; i < cumulativeProbability.size(); ++i) {
    for (unsigned j = 0; j < i; ++j) {
      cumulativeProbability[i] += probability[j];
    }
  }

  // испытания
  const int numberOfTests = 1000000;
  for (int i = numberOfTests; i > 0; --i) {
    double point = rnd() * cumulativeProbability.back(); // генерация числа от 0.0 до
                                                         // cumulativeProbability.back()
    int    indexRandomVariable = 0;
    for (; point > cumulativeProbability[indexRandomVariable]; ++indexRandomVariable);
    ++numberOfHits[indexRandomVariable - 1];
  }

  vector<double> experimentalProbability(len);
  for (unsigned i = 0; i < experimentalProbability.size(); ++i) {
    experimentalProbability[i] = numberOfHits.at(i) / (double)numberOfTests;
  }

  // результаты
  ofstream in;
  if (argc != 2) {
    in.open("graph");
  } else {
    in.open(argv[1]);
  }

  for (unsigned i = 0; i < experimentalProbability.size(); ++i) {
    in << probability.at(i) << " " << experimentalProbability.at(i) << endl;
  }

  in.close();
  return 0;
}

/*
  Биномиальное распределение
  Вход: индекс элемента распределения, общее число элементов
  Выход: вероятность элемента
*/
long double binomialDistribution(const int k, const int n, const double p) {
  double const q   = 1.0 - p;
  long double  res = 1.0;
  int          i   = 1;

  // Оууу, попробуй разберись
  for (; i <= k; ++i) {
      res *= p * (n - k + i) / i;
      if (i <= n-k) res *= q;
  }
  for (; i <= n-k; ++i) res *= q;

  return res;
}

/*
  Генератор случайных чисел [0.0,... 1.0]
*/
inline double rnd() {
  return rand() / (double)RAND_MAX; // от 0.0 до 1.0
}
